Bounds for the phonon-roton dispersion 
in superfluid 4 He 



J. Boronat a , J. Casulleras a , F. Dalfovo b , S. Moroni c , and S. Stringari 6 
a Departament de Fisica i Enginyeria Nuclear, Campus Nord B4-B5, 
Universitat Politecnica de Catalunya, E-08028 Barcelona, Spain 
b INFM, Dipartimento di Fisica, Universitd di Trento 1-38050 Povo, Italy 
c Istituto Nazionale di Fisica della Materia, Laboratorio FORUM, 
Scuola Normale Superiore, Piazza dei Cavalieri 7, 1-56126 Pisa 
(February 28, 1995) 

Abstract 

The sum rule approach is used to derive upper bounds for the dispersion law 
uJo(q) of the elementary excitations of a Bose superfluid. Bounds are explicitly 
calculated for the phonon-roton dispersion in superfluid 4 He, both at equilib- 
rium = 0.02186 A~ 3 ) and close to freezing (p = 0.02622 A~ 3 ). The bound 
wq(q) — 25(g) | x(q) where S(q) and x{l) are the static structure factor 
and density response respectively, is calculated microscopically for several val- 
ues of the wavevector q. The results provide a significant improvement with 
respect to the Feynman approximation ujfW) = q 2 (2mS(q))~ 1 . A further, 
stronger bound, requiring the additional knowledge of the current correlation 
function is also investigated. New results for the current correlation function 

are presented. 
67.40.-w, 67.40.Db 
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I. INTRODUCTION 



The microscopic investigation of the dynamic behavior of superfluid 4 He has been the 
object of extensive theoretical work in the past starting from the pioneering works by Bijl and 
Feynman [|l|] (see, for instance, Ref. [Q] for an up-to-date review). Recent approaches, based 



have provided accurate predictions for the dispersion law of this strongly interacting Bose 
system. 

The purpose of this paper is to show that in a Bose superfluid the knowledge of relevant 
static properties of the system can be used to derive useful upper bounds for the excitation 
spectrum, employing a sum rule approach (for a recent discussion on sum rules in Bose 
superfluids see for example Ref. ||). We will show that a key role in this context is played by 
the static density response for which Diffusion Monte Carlo (DMC) calculations have recently 
become available ||. Another relevant quantity in this context is the kinetic structure 
function, for which new DMC results will be presented. 

Explicit results for various bounds at the equilibrium density, p = 0.02186 A -3 , as well 
as close to freezing, p = 0.02622 A~ 3 , will be given in the first part of the work. In the 
second part we will discuss in detail the behaviour of the static response function, extending 
the analysis of Ref. |7| to lower g's and high pressure. 



The most famous estimate of the dispersion law in a Bose superfluid was proposed many 
years ago by Bijl and Feynman [[[]]. The resulting dispersion can be written in the form 



on perturbation theory with correlated functions || and on the use of shadow variables M , 



II. BOUNDS FOR THE PHONON-ROTON DISPERSION 



u F (q) 



mi(q) q 2 



(1) 



mo(q) 2mS(q) 



where 




(2) 



2 



are the fcth-moments of the dynamic structure function S(q,u) (H = 1 in this work). In 
deriving Eq. ([I]) one has evaluated the moment mi through the well known f-sum rule || 

Ml) = l(\p-*[H,pJi)) = N^- (3) 

holding for systems of particles interacting with velocity independent potentials. The brack- 
ets (...) denote ground state averages, while H is the iV-body Hamiltonian of the system 

N V7 2 N 

ff=-E^ + I>(r-«), (4) 

1=1 1<J 

and p q is the density fluctuation operator 

N 

p q = E^ q - ri • (5) 

i=l 

The moment mo has been instead expressed in terms of the static structure factor S(q) 
through the equation 

Ml) = (P-qPq) = NS(q) . (6) 

Both results (^) and (0) have been obtained using the completeness relation. The static 
structure factor S(q) is known with great accuracy both from Monte Carlo calculations and 
experimental data. In the present work we use the Diffusion Monte Carlo results shown in 
Fig. |. 

The Feynman energy (|l|), being based on the ratio of the two moments mi and m , 
provides, at zero temperature, a rigorous upper bound to the energy u (q) of the lowest 
state excited by the density operator pt, 

w (g) < uf(q) ■ (7) 

In the following we will identify the energy Uo(q) with the one of the elementary excitations 
of the system, i.e., the phonon-roton spectrum (in a Bose superfluid this identification is 
exact apart from decay processes of the elementary mode into two or more excitations ||). 
The bound (fj) reproduces exactly the phonon dispersion at small q: 
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u {q) = cq (8) 

where c is the sound velocity. This follows from the low q behavior of the static structure 
factor 

and is the consequence of the fact that in the macroscopic regime both the moments mo and 
mi are exhausted by the phonon mode. 

At higher wave vectors the Feynman bound instead overestimates significantly the ex- 
perimental dispersion law (see Fig. ||]). This is due to the occurrence of multipair excitations, 
whose sthength distribution, at energy higher than uj (q), turns out to be particularly im- 
portant in the determination of the energy weighted moment. For example the experimental 



data of Ref. |TT| at s.v.p. indicate that the roton exhausts only 1/3 of the energy weighted 
sum rule, the remaining part being associated with high energy multipair excitations. 

The idea to go beyond the Feynman approximation employing a sum rule approach was 
first developed many years ago by Feenberg |T^] with the help of the moments m 2 and 
777.3. Here we show that better bounds can be calculated using the inverse energy weighted 
moment 

m^q) = f du ^1 (10) 

This moment is an ideal quantity in order to investigate the collective properties of a Bose 
superfluid. In fact the factor 1/uj quenches significantly the high frequency tail of S(q,u) 
where multipair excitations are important. Furthermore the absence of single particle exci- 
tations in the low energy part of the spectrum, typical feature of a Bose system, makes the 
integral ( |T0"D particularly sensitive to the contribution of the collective mode. For the same 
reason the experimental determination of m_i, through a direct integration of the dynamic 
structure factor measured by neutron scattering, turns out to be more accurate than the 
one of any other moment [JTTJ . 
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The inverse energy weighted sum rule fllOD is directly related to the static density response 
of the system through the equation 

X (q) = -2m_i(g) (11) 

The static response fixes the linear changes in the density induced by an external static field 
interacting with the system with a potential of the form H ext = Ap q , coupled to the density 
fluctuation operator p q . At small q it yields the compressibility of the system 

X(0) = ~ (12) 
mc 2 

while at larger q is characterized by the occurrence of a pronounced peak (see discussion 
in Sect. HI). The static response x(<?) has been recently calculated in superfluid 4 He using 
Diffusion Monte Carlo techniques ||. These calculations reproduce the experimental data 
of x a ^ s.v.p. with good accuracy. 

The knowledge of the static response, together with the one of the static structure factor 
can be used to calculate a new upper bound for the dispersion law using the ratio 

mo(9) 25(g) 
^0-1(9) = = I 7 \ I 13 

between the non energy weighted and the inverse energy weighted sum rules. Since at zero 
temperature S(q,ui) = for uj < 0, the following inequality rigorously holds: 

w (g) < ^0-1(9) < wf(q) ( 14 ) 

At small q also the bound ( pT3[) approaches the phonon dispersion law as one can immediately 
see using results (§) and ( jl2[) . In Fig. || the new bound is reported for several values of q. 
The improvement with respect to the Feynman bound is significant both in the roton and 
in the maxon region. As already anticipated this improvement is the consequence of the 
fact that the moments mo and m_i, entering Eq. ([13]), are much less affected by multipair 
excitations with respect to the moment mi. This is also true at high pressure, as shown in 
Fig- |. 



A further improvement of the bound ([13]) can be obtained with the help of the energy 
weighted moments mi and m 2 . In fact one can derive the following inequality for the 
excitation energy uJo(q): 



^o(g) < 2 



^o-i - e) 2 + 4cj -iA 



(15) 



This bound is stronger than the bound (M). It involves the knowledge of the variance 

mi{q) m (q) 



A(g) 



m (q) m_i(g) 



(16) 



and of the energy 



e(q) = A- L (q) 



m 2 (q) 


( m (q) \ 


m (q) 


\m-i(q)J 



mi(q) 



(17) 



m_i(g) 

The latter depends not only on the moments m_i, mo and mi, already discussed above, 
but also on the moment m 2 . The quantity A vanishes when g — ^ since the two energies 
^o-i (<?) an d Uf(q) coincide in this limit as already pointed out before. Viceversa, the energy 
e, which represents an average energy of multipair excitations, is expected to depend less 
critically on q. Inequality (p~5| ) can be derived by using the fact that S(q, uj) vanishes for uj 
less than u>o(q) in bulk liquid 4 He; thus the following inequality holds for any positive 7: 

fduS(q,u)(l + 7cu) 2 



(18) 



J dujS(q, uj)uj~ 1 {1 + 7a;) 2 
The same inequality can be written in terms of the moments m^, and the value of 7 can 
be chosen in such a way to minimize the right hand side. After some algebra one obtains 
inequality (|T5|). A similar procedure can be used to derive upper and lower bounds to the 
static response function (see Ref. 0), and will be employed in Sect. PH . 

The moment m 2 was first explored by Feenberg |12] and turns out to be proportional to 
the current correlation function. In fact one has 



m 2 (q) = q 2 (J\jz q ) 



(19) 



where J zq is the z-component of the current density operator, and q is taken in the z- 



direction. Using the definition of the kinetic structure function |12 
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D (<?) = [ N i i ]\ f dr i • • • dr N [cos(q ■ r 12 ) - l](q ■ V^ )(q ■ V 2 ^ ) (20) 
the following expression for the moment rri2(q) holds: 

^ = f ( 2- S(q)) + lD (q) (21) 
mi{q) 2m m 

The kinetic structure function and, hence, the moment m 2 , has been directly calculated 
using a Diffusion Monte Carlo alghoritm. The results are shown in Figs. [| and |5J The 
structure of D(q) is almost the same at the two densities here considered; the curve under 
pressure is shifted upwards. This is just what one expects by looking at the large q behaviour 
of definition (|20|) , which yields D(q) — > (2iii/3)(Ek), for q 3> 27rp 1//3 . Our values of the mean 
kinetic energy (E K ) are 14.32(5) K and 19.57(5) K at p = 0.02186 A" 3 and p = 0.02622 
A~ 3 respectively. They give an asymptotic shift of 0.3 A~ 2 , in agreement with the data 
plotted in Fig. |j. The curve for D(q) at equilibrium density is also similar to the one used 
in Ref. J7j; in that case, the quantity D(q) was obtained by Fourier transforming the results 
of Path Integral Monte Carlo calculations [13] of the current correlation function in r-space. 

The microscopic results for the moments mo, m_i and m 2 allow one to calculate the 
bound (|I~5D; the results are reported in Figs. |2] and ||. All the moments entering this analysis 
have been calculated employing the Aziz potential HFDHE2 [|l6j. One notices a systematic 
improvement with respect to the bound (|13|) in the whole range of wavelength from maxons 
to rotons. At equilibrium density the bounds (HD, ( pT3|) and (|15|) yield the roton minimum at 
about 17.5K, 11.8K and 10. 8K, respectively, to be compared with the experimental value 
8.6K. The error bars in the figures are due to statistical errors in the calculation of m_i and 
m 2 . As concerns the pressure dependence, we note that the roton minimum shifts slightly 
to higher wave vectors by increasing pressure, in agreement with the experimental trend 



T4| . Also the roton gap exhibits the correct trend, being smaller at p = 0.02622. However, 



the statistical error on m_i prevents an accurate comparison with the experimental shift; in 
fact, the experimental roton gap decreases by only 1.3 K in the same range of pressure. 
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III. BOUNDS FOR THE STATIC RESPONSE FUNCTION 



So far we have used theoretical data for the moments m_i, m , mi and m 2 in order 
to evaluate rigorous upper bounds for the phonon-roton dispersion. Here we apply the 
formalism of Ref. [0] to evaluate upper and lower bounds to the moment m_i(g), and hence 
to the static response function x{q)i using mo, mi, m 2 , m 3 , as well as the experimental 
phonon-roton dispersion. This procedure will provide a check of consistency between the 
available theoretical calculations for the moments rrik{q), extending the analysis of Ref. 0] 
to lower g's and to high pressure. 

One can easily derive lower bounds to m_i(g) starting from the inequality 

r dLO^^-il + aoo + /V) 2 > (22) 
Jo us 

holding for any real a and (3. The same inequality can be written as a lower bound to m_i(q). 
Minimization with respect to the parameters a and (3 provide the bounds. In particular, by 
minimizing with respect to a with [3 = one gets the Feynman approximation to m_i(q): 

m_i(?) > m F x (g) = 2Nmq~ 2 S 2 (q) . (23) 

Minimizing with respect to both a and [3 one gets a stronger lower bound: 

--^-T^mk) • (24) 

where 

= mM rnM (25) 
mi{q) m [q) 

and 



e(q) = A' 1 



m 3 (q) . ( m x {q) 2 o m 2 (q) 



(26) 



m x (q) mo(q) m (q)_ 

One notes that the quantities A and e have the same form of A and e in Eqs. ( |16D and (|17|), 
but with the index of the k- moments scaled by 1. Inequality (E3) requires the knowledge 



S 



of the cubic energy weighted moment 7713(g), which can be calculated through the Puff sum 



rule [17 



m 3 (g) = N 



2m rnZ 



11 
777/ 



2m 2 



J dv g(r)(l — cos(q-r))(q- V) 2 V{r) 



(27) 



where V(r) and g(r) are the interatomic potential and the radial distribution function, 
respectively. 

In a similar way, one can derive upper bounds 0. One finds 



m_ 1 (g)<^M = iV5(gV - 1 (g) 



as well stronger upper bound: 



m-i(q) < 



m (q) 
"oil) 



m (q) ( mx{q) 



' m 2 {q) 
mi(g) 



^o(?) 



(28) 



(29) 



™a(g) \mo(q) 

The results for the above lower and upper bounds for m_ 1 (g) are shown in Figs. |6| and 
|7| at equilibrium density and close to freezing, respectively. To evaluate the bounds we have 
used the same moments as in Sect. fTI] and the experimental phonon-roton dispersion 
for u; (g). Dashed lines correspond to the weakest bounds (p3[) and (p8|) , while solid lines 
correspond to the bounds (|24) and (p9|). The latter account for the effect of multiphonon 
excitations through the inclusion of higher /c-moments. This explains why the allowed area 
for m_i, i.e. between lower and upper bounds, is significantly reduced passing from dashed 



to solid lines. Indeed the bounds (|24]) and ( p9|) represent a quite stringent test of consistency 
between independent calculations and measurements of fc-moments of the dynamic structure 
function S(q,u). The available experimental data of m_x at equilibrium are consistents 
with the bounds. The same is true for the Diffusion Monte Carlo data at equilibrium 
and freezing pressure. We note that the new Monte Carlo data for 5(g) and -D(g) provides 
accurate bounds even at relatively small g's, i.e, in the maxon region 0.5 to 1.5 A" 1 . 



IV. CONCLUSION 



In the first part of this work we have discussed new upper bounds for the excitation 
spectrum in superfluid 4 He. The method makes use of basic static properties of the system: 
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the static structure factor, the static response and the current correlation function. These 
quantities are now available in microscopic ab initio calculations with good accuracy. In 
particular we have used recent Diffusion Monte Carlo data for the static response || and we 
present new results for the kinetic structure function. The upper bounds for the phonon- 
roton dispersion turn out to be rather close to the experimental values and can be calculated 
at any pressure. 

In the second part we have evaluated upper and lower bounds to the static response 
function using a method proposed by two of us in Ref. 0. The new data for the current 
correlation function allows one to extend the analysis of Ref. to lower values of q and to 
high pressure. The main result is a general consistency between the independent evaluations 
of the several fc-moments involved in the analysis, and hence of the quantities S(q), D(q), 
x{q) in a wide range of q. 
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FIGURES 

FIG. 1. Static structure factor S(q) at equilibrium density (solid line) and at p = 0.02622 A -3 
(dashed line) 

FIG. 2. Phonon-roton spectrum at the equilibrium density, p = 0.02186 A -3 . Solid line: 
experiments [jOJ; dashed line: Feynman approximation (|l]); empty circles: upper bound Wo-l 
defined in Eq. (|i~3|); solid circles: upper bound (|l5|). 



FIG. 3. Same as in Fig. || but at freezing pressure (p = 0.02622 A" 3 ). The solid line 



corresponds to a smooth interpolation between experimental data on the phonon dispersion [13], 
up to 1 A -1 , and recent data on the roton minimum [14|. 



FIG. 4. Kinetic structure function D(q) at equilibrium density (empty circles) and close to 
freezing (solid circles). 

FIG. 5. Ratio mi(q)/m\(q) at equilibrium density (empty circles) and close to freezing (solid 
circles) 

FIG. 6. Inverse energy weighted moment m-\{q) at equilibrium density. Empty circles: 
experiments [|ll|]; solid circles with error bars: Diffusion Monte Carlo calculations ||; dashed lines: 
upper and lower bounds fl28|) and (|23l); solid lines: upper and lower bounds (|29|) and fl2~3|). 



FIG. 7. Inverse energy weighted moment m-\{q) at density p = 0.02622 A 3 . Points with 
error bars: Diffusion Monte Carlo calculations ||; dashed lines: upper and lower bounds (pq) and 



|23|) ; solid lines: upper and lower bounds fl29[) and (|24|) . 
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